Accessing & Mapping OSM in R

Setup

Before we begin our exercise, we need to install all the needed packages - we only need to use install.packages() once for each. However EVERY time you start using R you will need to run library() to load the installed package.

In this case we will use if(!require("osmdata")) install.packages("osmdata") to check if the package or library called osmdata is installed, if it is, we load it, and if it’s not it will be installed for us.

In this exercise we are going to leverage the package osmdata to directly query and download features from OSM. This package will directly access the Overpass Turbo API to create and make queries directly from the R environment. Nevertheless, the use of the overpass-turbo.eu can be useful when we are not sure what we are looking for or when we have some difficulty in building the query.

See here for some basics on OverpassTurbo and a formal guide on queries.

Building a OSM query

In order to structure all the data in the world, OSM created a hierarchical data structure. In this case we can think of the world being broken down into spatial features (e.g. points, lines, polygons) each with a tag made up of a Key (e.g. amenity, building, highway etc) and Value (e.g. coffee shop, primary road).

We can look the tagging conversion on OSM’s map features directory.

We can list all the the available OSM features using available_features() to list the available feature keys. In this case - to limit the length of the output we will use head() so we only view the first 5 features listed.

# list top five
head(available_features())

## [1] "4wd_only"  "abandoned" "abutters"  "access"    "addr"      "addr:city"

# list all
# available_features()

Now we can see what the different types of buildings are based on their tag value. Here again we are using head() so that we only see the first five results.

# list first 5 building types
head(available_tags("building"))

## [1] "apartments" "bakehouse"  "barn"       "barracks"   "beach_hut" 
## [6] "bridge"

# list all building types
# available_tags("building")

————————————————————–

Challenge

Use the available_features() and available_tags("a_feature_name") to find the feature and tag name for stores that sell alcohol.

# enter your answer here

————————————————————–

Your First Query - Exploring Jamaica

In this example we will be exploring the area around you in Jamaica. First it helps to go to OSM https://www.openstreetmap.org/ and figure how features and tags can be found.

To find feature tag pairs zoom to a business of interest, right-click > Show Address > click on link to results > top row is [feature, tag] in this case [amenity, bank].

Finding a feature and tag on OSM

Ok let’s start by trying to find and download data on restaurants in Montego Bay. To do this we need to create a bounding box using getbb(). Here make sure to be explicit, let’s search for "Montego Bay, Jamaica" to avoid other locations called Montego Bay. This then gets passed to the opq() function to convert that into an overpass query based on location, then we pass an additional argument to the query using add_osm_feature(), in particular we are going to look for, in this case a feature key "amenity" and a tag value "restaurant".

query = getbb(place_name="Montego Bay, Jamaica") %>%
      opq() %>%
       add_osm_feature(key="amenity", value="restaurant")

We can then look at how the query is structured before it is passed to overpass turbo.

query

## $bbox
## [1] "18.4630457,-77.9274209,18.4790101,-77.9119878"
## 
## $prefix
## [1] "[out:xml][timeout:25];\n(\n"
## 
## $suffix
## [1] ");\n(._;>;);\nout body;"
## 
## $features
## [1] " [\"amenity\"=\"restaurant\"]"
## 
## attr(,"class")
## [1] "list"           "overpass_query"
## attr(,"nodes_only")
## [1] FALSE

Although it looks complex, we can see it has a few important components, one called bbox which is the bounding box (extent of) for Montego Bay, it as features with key "amenity" and tag value "restaurant".

We can then send our query and retrieve the objects as “simple feature” sf objects in R. You can think of sf as being a component of R developed to handle spatial vector features (points, lines, polygons) and their attribute tables. For more on sf see here.

Run the query, depending on the size of the query and your internet speeds this might take a bit of time.

restaurant = osmdata_sf(q=query)
restaurant

## Object of class 'osmdata' with:
##                  $bbox : 18.4630457,-77.9274209,18.4790101,-77.9119878
##         $overpass_call : The call submitted to the overpass API
##                  $meta : metadata including timestamp and version numbers
##            $osm_points : 'sf' Simple Features Collection with 66 points
##             $osm_lines : NULL
##          $osm_polygons : 'sf' Simple Features Collection with 15 polygons
##        $osm_multilines : NULL
##     $osm_multipolygons : NULL

restaurant is now an object with a few different components restaurant$bbox will give us back the bounding box, restaurant$osm_points contains 66 point features (i.e. the restaurants) and restaurant$osm_polygons has 15 polygons (?the restraurants in polygon?).

To help us understand what we got back from overpass turbo, let’s put them on a map. To do this we will use the mapview function from the excellent mapview package. mapview will automatically generate an interactive webmap based on leaflet.

Plotting on a Map

Let’s create an interactive map by passing the restaurant$osm_points to mapview. Note we are turning off the legend by setting it to F, if you want it on, set legend=T.

mapview(restaurant$osm_points, legend=F)

mapview makes it extremely easy to combine multiple sets of data to a map. In this case we will just + two mapview objects of interests.

In this case let’s see what the difference between the restaurant$osm_polygons and restaurant$osm_points are. Zoom in an see if they are different.

mapview(restaurant$osm_polygons, legend=F) + 
  mapview(restaurant$osm_points, legend=F)

Looks like the points are just the exterior of our polygons. Maybe we just plot the polygons instead of the points.

Before we go on, try hovering over one of the points or polygons. Note that there is a popup window, with info that isn’t very informative

Let’s clean this up. To do this we need to track down the restaurant names somewhere in our data. Let’s see what other data is in the restaurant$osm_polygons data frame.

head(restaurant$osm_polygons)

## Simple feature collection with 6 features and 16 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -77.9245 ymin: 18.47092 xmax: -77.9226 ymax: 18.47656
## Geodetic CRS:  WGS 84
##              osm_id                 name addr.city addr.housenumber
## 368830821 368830821         Island Grill      <NA>             <NA>
## 371393488 371393488            Pizza Hut  Kingston             <NA>
## 371444457 371444457          Lee's Pizza      <NA>             <NA>
## 371444946 371444946 R&D Restaurant & Bar      <NA>             <NA>
## 371544807 371544807 Yash Bowl Restaurant      <NA>             <NA>
## 371545035 371545035   Chinese Restaurant      <NA>             <NA>
##                      addr.street    amenity building cuisine delivery level
## 368830821                   <NA> restaurant      yes    <NA>     <NA>  <NA>
## 371393488 Howard Cooke Boulevard restaurant      yes   pizza      yes     0
## 371444457                   <NA> restaurant      yes   pizza     <NA>  <NA>
## 371444946                   <NA> restaurant      yes    <NA>     <NA>  <NA>
## 371544807                   <NA> restaurant      yes    <NA>     <NA>  <NA>
## 371545035                   <NA> restaurant      yes chinese     <NA>  <NA>
##             name.it                                       opening_hours
## 368830821      <NA>                                                <NA>
## 371393488      <NA> Mo-Th 10:00-23:00; Fr-Sa 10:00-23:59;Su 08:00-23:00
## 371444457      <NA>                                                <NA>
## 371444946      <NA>                                                <NA>
## 371544807 Yash Bowl                                                <NA>
## 371545035   Chinese                                                <NA>
##                        operator           phone                source takeaway
## 368830821                  <NA>            <NA> KG Ground Survey 2015     <NA>
## 371393488 Restaurant of Jamaica 1(876) 971-5379 KG Ground Survey 2015      yes
## 371444457                  <NA>            <NA> KG Ground Survey 2015     <NA>
## 371444946                  <NA>            <NA> KG Ground Survey 2015     <NA>
## 371544807                  <NA>            <NA> KG Ground Survey 2015     <NA>
## 371545035                  <NA>            <NA> KG Ground Survey 2015     <NA>
##                                 geometry
## 368830821 POLYGON ((-77.92386 18.4747...
## 371393488 POLYGON ((-77.92444 18.4757...
## 371444457 POLYGON ((-77.92368 18.4765...
## 371444946 POLYGON ((-77.9227 18.47478...
## 371544807 POLYGON ((-77.92314 18.4710...
## 371545035 POLYGON ((-77.92402 18.4725...

Looks like there is a column called name which holds the name of each restaurant polygon. Let’s see if we can figure out how to access it. In R $ is often used to access data by column name. Let’s try restaurant$osm_polygons$name .

restaurant$osm_polygons$name

##  [1] "Island Grill"                   "Pizza Hut"                     
##  [3] "Lee's Pizza"                    "R&D Restaurant & Bar"          
##  [5] "Yash Bowl Restaurant"           "Chinese Restaurant"            
##  [7] "Union Street Restaurant & Cafe" "Phoenix Restaurant"            
##  [9] "Captain's Bakery & Grill"       "Stanley's Patties"             
## [11] "Nyam & Jam Restaurant"          "Hunger Buster"                 
## [13] "Homestyle Restaurant"           NA                              
## [15] "Pier 1"

Perfect. Got it.

Ok now let’s pass a few more arguments to the mapview function to pretty things up. We can set the on-click data with label the color of the polygon outline with color and the polygon fill color with col.regions.

mapview(restaurant$osm_polygons, 
        legend=F, 
        label = restaurant$osm_polygons$name, 
        color = 'red', 
        col.regions='red')  

————————————————————–

Challenge

Create a map containing the polygons for both restaurants and pharmacies in Kingston.

Remember in coding, its always recommend to copy and paste a working example.

# insert your code here

————————————————————–

Additional Tricks

We can also use grep to search the available features. In this case I will find any features that contain the word 'building’.

# search for the following pattern
search_string = 'building'

# search function (don't edit)
available_features()[grep(search_string, available_features(),ignore.case=T)]

## [1] "building"             "building:colour"      "building:fireproof"  
## [4] "building:flats"       "building:levels"      "building:material"   
## [7] "building:min_level"   "building:part"        "building:soft_storey"

————————————————————–

Other resources

#https://dominicroye.github.io/en/2018/accessing-openstreetmap-data-with-r/

# https://towardsdatascience.com/calculating-building-density-in-r-with-osm-data-e9d85c701e19

#https://dominicroye.github.io/en/2018/accessing-openstreetmap-data-with-r/

#https://docs.ropensci.org/osmdata/